function [z, df, d2f] = FUNC1(a,b,c,d,e,f)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This version is checked on Oct 30, 2023
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% a omega0_t
% b RL_t+1
% c H_t
% d RK_t+1
% e sigma0
% f sigma1

%checked
fun = @(x,y)exp(x+y).*exp(-(x+0.5.*e.^2).^2/(2.*e.^2))./(e.*sqrt(2.*pi)).*exp(-(y+0.5.*f.^2).^2./(2.*f.^2))./(f.*sqrt(2.*pi));
ymin = @(x)log(b.*c./(d.*(1+c)))-x;

z = integral2(fun,a,100,ymin,100);

%cheacked
funa =@(s)exp(a+s).*exp(-(a-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi)).*exp(-(s-(-0.5.*f.^2)).^2./(2.*f^2))./(f.*sqrt(2.*pi));
da=-integral(funa,log(b*c/(d*(1+c)))-a,100);

%cheacked
funb = @(x)c./(d.*(1+c)).*exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi));
db=-integral(funb,a,100);

% checked
func = @(x)b./(d.*(1+c).^2).*exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi));
dc=-integral(func,a,100);

% checked
fund = @(x)(c.*b)./(d.^2.*(1+c)).*exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))/(e.*sqrt(2.*pi));
dd=integral(fund,a,100);

% checked
fune = @(x,y)exp(x+y).*(-exp(-(x+0.5.*e.^2).^2/(2.*e.^2))./(e.^2.*sqrt(2.*pi))+exp(-(x+0.5.*e.^2).^2./(2.*e.^2))./(sqrt(2.*pi)).*(-((x+0.5.*e.^2).*e.^2-(x+0.5.*e.^2).^2)./e.^4)).*exp(-(y+0.5.*f.^2).^2./(2.*f.^2))./(f.*sqrt(2.*pi));
de = integral2(fune,a,100,ymin,100);

%checked
fung = @(x,y)exp(x+y).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2))./(e.*sqrt(2.*pi)).*(-exp(-(y+0.5.*f.^2).^2./(2.*f.^2))./(f.^2.*sqrt(2.*pi))+exp(-(y+0.5.*f.^2).^2./(2.*f.^2))./(sqrt(2.*pi)).*(-((y+0.5.*f.^2).*f.^2-(y+0.5.*f.^2).^2)./f.^4));
dg = integral2(fung,a,100,ymin,100);

df=[da db dc dd de dg];

% second order derivatives
% a omega0_t
% b RL_t+1
% c H_t
% d RK_t+1
% e sigma0
% f sigma1
    % define the probablity density function
    function [f0] = desity(x,sigma_under)
      f0 = 1/(sigma_under*sqrt(2*pi))*exp(-(x+0.5*sigma_under^2)^2/(2*sigma_under^2));
    end
    
    % call the probability density function
    f0_omega0bar = desity(a,e);
    f1_omega0bar = desity(log(b*c/(d*(1+c)))-a,f);
    
    % first order derivative of probability density function with respect
    % to sigma
    function [f0prime] = desityprime(x,sigma_under)
      f0prime = -1/(sigma_under^2*sqrt(2*pi))*exp(-(x+0.5*sigma_under^2)^2/(2*sigma_under^2))+ 1/(sqrt(2*pi))*exp(-(x+0.5*sigma_under^2)^2/(2*sigma_under^2))*(-(((x+0.5*sigma_under^2)*sigma_under^2 -(x+0.5*sigma_under^2)^2)/(sigma_under^4)));
    end
    
    % Second order derivative of probability density function
    function [f0prime2] = desityprime2(x,sigma_under)
      f0prime2 = 2/(sigma_under.^3.*sqrt(2.*pi)).*exp(-(x+0.5.*sigma_under.^2).^2./(2.*sigma_under.^2)) ...
               + 1/(sqrt(2*pi)).*exp(-(x+0.5.*sigma_under.^2).^2./(2.*sigma_under.^2)).*((((x+0.5.*sigma_under.^2).*sigma_under.^2 -(x+0.5.*sigma_under.^2).^2)/(sigma_under.^5)))...
               + 1/(sqrt(2*pi)).*exp(-(x+0.5.*sigma_under.^2).^2./(2.*sigma_under.^2)).*(-((2.*x.*sigma_under.^2+2.*sigma_under.^4-6.*sigma_under.^2.*(x+0.5.*sigma_under.^2)+4.*(x+0.5.*sigma_under.^2).^2)./(sigma_under.^5)))...
               + 1/(sqrt(2*pi)).*exp(-(x+0.5.*sigma_under.^2).^2./(2.*sigma_under.^2)).*(-(((x+0.5.*sigma_under.^2).*sigma_under.^2 -(x+0.5.*sigma_under.^2).^2)./(sigma_under.^4))).*(-(((x+0.5*sigma_under.^2).*sigma_under.^2 -(x+0.5*sigma_under.^2).^2)/(sigma_under.^3)));
    end
    f0prime_omega0bar = exp(-(a-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi)).*(-(a-(-0.5.*e.^2))./(e.^2));

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For da functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %checked
    funaa = @(x)exp(a+x).*(f0_omega0bar+f0prime_omega0bar).*exp(-(x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi));
    daa1=integral(funaa,log(b*c/(d*(1+c)))-a,100);
    % get daa
    daa = -(b*c/(d*(1+c))).*f0_omega0bar.*f1_omega0bar-daa1;
    
    % get dab %checked
    dab = (b*c/(d*(1+c))).*f0_omega0bar.*f1_omega0bar./b;
    
    % get dac %checked
    dac = (b*c/(d*(1+c))).*f0_omega0bar.*f1_omega0bar./(c.*(1+c));
    
    % get dad %checked
    dad = (b*c/(d*(1+c))).*f0_omega0bar.*f1_omega0bar./(-d);
    
    % get dae %checked
    f0prime_omega0bar_sigma0 = desityprime(a,e);
    funae = @(x)exp(a+x).*f0prime_omega0bar_sigma0.*exp(-(x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi));
    dae = -integral(funae,log(b*c/(d*(1+c)))-a,100);
    
    % get dag %checked
    funag = @(x)exp(a+x).*f0_omega0bar.*(-1./(f.^2.*sqrt(2.*pi)).*exp(-(x+0.5.*f.^2).^2./(2.*f.^2))+ 1./(sqrt(2.*pi)).*exp(-(x+0.5.*f.^2).^2./(2.*f.^2)).*(-(((x+0.5.*f.^2).*f.^2 -(x+0.5.*f.^2).^2)/(f.^4))));
    dag = -integral(funag,log(b*c/(d*(1+c)))-a,100);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For db functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % get dba 
    dba = dab; %checked
    
    % get dbb % checked
    funbb = @(x)(c./(d.*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./b;
    dbb = -integral(funbb,a,100);
    
    % get dbc  %checked
    funbc1 = @(x)(c./(d.*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2))./(f.^2)))./(c.*(1+c));
    funbc2 = @(x)(1./(d.*(1+c).^2)).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)));
    
    dbc = -integral(funbc1,a,100)-integral(funbc2,a,100); 
    
    % get dbd % checked
    funbd1 = @(x)(c./(d.*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(-d);
    funbd2 = @(x)(c./(d.^2.*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)));
    
    dbd = -integral(funbd1,a,100)+integral(funbd2,a,100);
       
    % get dbe % checked
    funbe = @(x)(c./(d.*(1+c))).*exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-1./(e.^2.*sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2))+ 1./(sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2)).*(-(((x+0.5.*e.^2).*e.^2 -(x+0.5.*e.^2).^2)./(e.^4))));
    dbe = -integral(funbe,a,100);
    
    % get dbf %checked
    funbg = @(x)(c./(d.*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi))).*(-1/(f.^2.*sqrt(2*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2))+ 1./(sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1+c)))-x+0.5.*f.^2).^2./(2.*f.^2)).*(-(((log(b.*c./(d.*(1+c)))-x+0.5.*f.^2).*f.^2 -(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2)./(f.^4))));
    dbg = -integral(funbg,a,100);
      
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For dc functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    
    % get dca %checked
    dca = dac;
    
    % get dcb %checked
    dcb = dbc;
    
    % get dcc; %checked, found a mistake on oct 30, 2023
    funcc1 = @(x)(b/(d*(1+c)^2)).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(c.*(1+c));
    funcc2 = @(x)(2*b/(d*(1+c)^3)).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)));

    dcc =  -integral(funcc1,a,100) + integral(funcc2,a,100);
    
    % get dcd %checked, found a mistake on oct 30, 2023
    funcd1 = @(x)(b./(d.*(1+c).^2)).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(-d);
    funcd2 = @(x)(b./(d.^2.*(1+c).^2)).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)));

    dcd = -integral(funcd1,a,100) + integral(funcd2,a,100);
    
    % get dce, %checked, found a mistake on oct 30, 2023
    funce = @(x)(b./(d.*(1+c).^2)).*exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-1/(e.^2.*sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2))+ 1./(sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2)).*(-(((x+0.5.*e.^2).*e.^2 -(x+0.5.*e.^2).^2)./(e.^4))));
    dce = -integral(funce,a,100);
    
    % get dcg, %checked, found a mistake on oct 30, 2023
    funcg = @(x)(b./(d.*(1+c).^2)).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(-1./(f.^2.*sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2))+ 1./(sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2)).*(-(((log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).*f.^2 -(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2)./(f.^4))));
    dcg = -integral(funcg,a,100);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For dd functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    
    % get dda %checked
    dda = dad;
    
    % get ddb %checked
    ddb = dbd;
    
    % get ddc; %checked
    ddc = dcd;
    
    % get ddd %checked
    fundd1 = @(x)(c*b/(d^2*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(-d);
    fundd2 = @(x)(2*c*b/(d^3*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)));

    ddd =  integral(fundd1,a,100) - integral(fundd2,a,100); 
    
    % get dde %checked
    funde = @(x)(b.*c./(d.^2.*(1+c))).*exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-1/(e.^2*sqrt(2*pi))*exp(-(x+0.5.*e.^2).^2./(2.*e.^2))+ 1/(sqrt(2.*pi)).*exp(-(x+0.5*e.^2).^2./(2.*e.^2)).*(-(((x+0.5.*e.^2).*e.^2 -(x+0.5.*e.^2).^2)./(e.^4))));
    dde = integral(funde,a,100);
    
    % get ddg % checked
    fundg = @(x)(b.*c./(d.^2.*(1+c))).*exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi)).*(-1./(f.^2.*sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2))+ 1/(sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2)).*(-(((log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).*f.^2 -(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2)./(f.^4))));
    ddg = integral(fundg,a,100);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For de functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    
    % get dea %checked
    dea = dae;
    
    % get deb %checked
    deb = dbe;
    
    % get dec %checked
    dec = dce;
    
    % get ded %checked
    ded = dde;
    
    % get dee; %checked
    % may be problematic
    funee = @(x,y)exp(x+y).*(desityprime2(x,e)).*exp(-(y+0.5.*f.^2).^2./(2.*f.^2))./(f.*sqrt(2.*pi));
    ymin = @(x)log(b.*c./(d.*(1+c)))-x;

    dee = integral2(funee,a,100,ymin,100);

    % get deg, %checked
    funeg = @(x,y)exp(x+y).*(-1/(e.^2.*sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2))+ 1./(sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2)).*(-(((x+0.5.*e.^2).*e.^2 -(x+0.5.*e.^2).^2)./(e.^4)))).*(-1./(f.^2.*sqrt(2.*pi)).*exp(-(y+0.5.*f.^2).^2./(2.*f.^2))+ 1/(sqrt(2.*pi)).*exp(-(y+0.5.*f.^2).^2./(2.*f.^2)).*(-(((y+0.5.*f.^2).*f.^2 -(y+0.5.*f.^2).^2)./(f.^4))));
    deg = integral2(funeg,a,100,ymin,100);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For dg functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    
    % get dga %checked
    dga = dag;
    
    % get dgb %checked
    dgb = dbg;
    
    % get dgc %checked
    dgc = dcg;
    
    % get dgd %checked
    dgd = ddg;
    
    % get dge %checked
    dge = deg;
    
    % get dgg %checked
    % may be problematic % becareful with this one
    fungg = @(x,y)exp(x+y).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2))./(e.*sqrt(2.*pi)).*(desityprime2(y,f));
    ymin = @(x)log(b.*c./(d.*(1+c)))-x;

    dgg = integral2(fungg,a,100,ymin,100);   
    
    
d2f = [daa dab dac dad dae dag; 
       dba dbb dbc dbd dbe dbg;
       dca dcb dcc dcd dce dcg;
       dda ddb ddc ddd dde ddg;
       dea deb dec ded dee deg;
       dga dgb dgc dgd dge dgg];

end



